Diagrammatic Quantum Monte Carlo Algorithm in Momentum Representation: 
Hess-Fairbank Effect and Mesoscopics in ID BEC with Attractive Interaction 
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A novel algorithm of Diagrammatic Quantum Monte Carlo in momentum representation is re- 
ported in details. New models can be studied with this algorithm. For Bose systems with attractive 
interaction, the algorithm is free of the well-known minus-sign problem, while in other models it is 
weaker than in real-space methods. 

Using this algorithm, we present the results of an exact numeric simulation of A'' one-dimensional 
bosons with attractive 5-functional interaction in a rotating ring. We prove that in the large-A limit 
the system can be described by conventional methods of weakly interacting gas, the dimensionless 
parameter of weak interaction being just 1/N . When the strength of interaction is less then a certain 
threshold value, the dependence of angular momentum on the rotation frequency features plateaus 
characteristic of the irrotational fluid (the Hess-Fairbank effect). 

PACS numbers: 03.75.Fi, 02.70.Ss, 68.65.-k 
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I. QMC ALGORITHM IN MOMENTUM 
REPRESENTATION AGAINST SIGN PROBLEM 

A. Background: sign problem 

The Quantum Monte Carlo (QMC) methods (varia- 
tional, determinant, trajectory etc.) prove their useful- 
ness for studying thermodynamics of diverse quantum 
systems. In QMC, the partition function Z = Tie~^^ 
is broken into series and summed using importance sam- 
pling. Each term is represented by unique set of inner 
parameters (MC configuration) and can be positive or 
negative depending on the actual MC configuration, the 
model Hamiltonian, and particle statistics. Alternating 
sign results in fluctuation of partition function and other 
quantities used in calculation. As the temperature is low- 
ered, the errors grow larger making the calculation slug- 
gish or even impossible. This is the so-called minus-sign 
problem, the inherent feature of most trajectory QMC 
methods. 

Some cluster methods are free of minus sign. The de- 
terminant QMC method |^] is free of minus sign, but un- 
fortunately, its application is restricted by fermionic and 
spin systems. Moreover, its running time as a function 
of cluster size L scales as while in trajectory methods 
it is linear by L. So for large clusters, the expected slow- 
ing down of determinant method is stronger than that of 
trajectory method caused by sign problem. The method 
of exact diagonalization of Hamiltonian |^ is more use- 
ful and applies to wider class of models. Nevertheless, 
the size of system is limited by Lmax — 10 12 as the 
amount of needed calculations grows exponentially with 
increasing L. 

Some sort of correction Q makes possible to determine 
the ground-state energy with good precision even if ap- 
proaching the temperature low enough is prevented by 
sign problem. For electrons on simple square (or cu- 
bic) lattice and hopping only to nearest neighbours, the 



particle-hole transormation helps to remove part of sign 
not linked with Fermi statistics. But generally, sign prob- 
lem is present in simulation. 

This Letter is organized as follows. In Section |, we 
present the trajectory algorithm in momentum repre- 
sentation developed in the framework of Diagrammatic 
QMC method weakening or removing minus-sign in many 
models. In Section we study the one-dimensional 
bosonic system with attractive interaction in a rotating 
ring taken by Ueda and Leggett as a model of irrota- 
tional fluid Q. This model can not be simulated with 
usual real-space methods, in part due to sign problem, 
while the algorithm described in this Letter is very effi- 
cient in this system allowing to simulate up to 100 and 
more particles. 



B. DQMC basics 

In this paragraph, we repeat shortly the DQMC basics 
to be used later [g. 

Decomposition. Diagrammatic QMC is based on the 
following decomposition of partition function into series 
of interaction representation: 

Z = ^ ^ / {~dTm-i) / (-dr„_2)... / (-dro) x 
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where the system Hamiltonian is H = Hq + V , energies 
Eq{...) are given by Hq which is chosen to be diagonal on 
occupation numbers {n}, Tm = tq-\- (3, {n^™^} = 

Each term in Eq. (j^) is represented by its own picture 
of particle trajectories in (x, r)-space, where r = /3, 
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f3 = l/fcsT. The weight of a given MC configuration 
writes as 



Wmc 



n (- Ar (... I F|...)exp (...)) 



Then mark the points of worldhne distortions, the so- 
called "kinks" 1^. The imaginary-time step Ar can be 
taken small enough (~ \0~^(3). 

Process of calculation. The importance sampling 
(Metropolis algorithm) consists in random transforma- 
tions of MC configurations obeying the following require- 
ments: 

• Full set. Two arbitrary MC configurations with 
nonzero weight can be transformed into each other 
with nonzero probability and in finite number of 
steps; 

• Balance. For each updating process transforming 
MC configurations A to B ("direct" process), we 
juxtapose the "inverse" process transforming MC 
configurations B exactly to A. Additionally, their 
frequencies must meet the balance equation 

I Wa I P^pt^"^ = \Wb\ P^pt^'^'i , 

with Wa and Wb weights of MC configurations A 
and B, P^ and P^ the frequencies to call direct and 
inverse processes, and p^^'^^ and p^^'^^ the proba- 
bilities to accept the respective update. The latter 
are usually determined from the relations 



^(acc) _ 



where a 



p. 



WaP^ 
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if a < 1, 
if a > 1, 



is the so-called "acceptance ra- 
tio". Therefore, to determine the probability of 
accepting the update, we must know parameters of 
respective inverse process. 

There remain some freedom of choosing the time r 
of newly-created (or shifted) kink. The most reasonable 
approach is to choose T with probability * 
using relation 



Zi{SE,Trr 



— '^min ^ f^(A-Z?, Tjjidx '^min) (2) 

(see Appendix ^) , according to the fact that the weight 
of new MC configuration is proportional to e~^^'^ . 



C. DQMC in momentum representation 

The algorithm described here is developed for systems 
with a Hubbard-like Hamiltonian in momentum repre- 
sentation 



E 

p,q,r,s 



(3) 




FIG. 1: Kinks generated by two-particle interaction: a, b - 
diagonal, c, d ~ non-diagonal; a, b in case of contact interac- 
tion Upqrs ~ UoSp+q=r+s Can be taken into account analyti- 
cally; a, c appear in simulation of Bose systems only. Here 
and troughout the paper, imaginary-time axis is plotted hor- 
izontally, momenta are placed vertically, occupation number 
indicated by line thickness. 
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FIG. 2: 
kinks. 



Sample Monte Carlo configuration with all types of 



Though typically, interaction term conserves full mo- 
mentum Upqrs = Upq5p+q^r+s, but it is uot mandatory 
for this algorithm. 

Kinks. The interaction term taken as aperturbation, 
generates four-ended kinks shown in Fig. The multi- 
plier entering the configuration weight due to each kink, 
is given by sum of respective Upqrs ^^npUqU^Us for all non- 
identical recombinations of momenta p,q,r,s. For Fermi 
systems, half of these terms get negative sign. This allows 

E (m^^'-i) 

to use the standard relation Sign (F) = (— 1) * 

linking fermionic sign of configuration and time wind- 

(r) 

ing numbers W^' of worldlines which holds in real-space 
trajectory methods. 

The example of MC configuration is sketched in Fig. |^. 
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(a) 




(c) 



(b) 



~| -•— () 
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FIG. 3: The updating processes for non-diagonal kinks: Cre- 
ation/Annihilation of Two Kinks (a), Expanding/Contracting 
of a Kink (b), Entangling of a Kink with a Worldline (c), and 
Shifting of a Kink Through Another (d). 
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FIG. 4: Creation of two "biased" kinks and annihilation of 
them in different way change full momentum of the system. 



of particles {'r]'^{a^ +a^), like in Worm algorithm^), or 

i 

full momentum {rj ^ a^flrOs, see Fig. with r] 

pqrs 
p+q=£r+s 

small enough. However, fictitious kinks are not needed 
when calculating energy levels characterized by their own 
value of full momentum. 



D. Applications of the method 



Sign problem. The sign problem in usual real-space 
QMC algorithm is caused by hopping term. In this al- 
gorithm, on the contrary, it is caused by sign of interac- 
tion. Thus, weakly-interacting systems can be simulated 
in momentum representation with much greater preci- 
sion. 

For attractive interaction, the algorithm is free of sign 
problem, as every kink enters MC weight with positive 
factor. The particle-hole transformation for electrons on 
a lattice can remove part of sign not linked with Fermi 
statistics, in case when their interaction is ^ Kj^^it'^ii- 

In simple but rather wide case of point interaction 
Upqrs — UoSp+q^r+si the diagonal part of interaction 
(kinks Fig. ^ (a,b)) is summed analytically and the MC 
configurations with single kink giving most part of sign 
problem, become impossible. Monte Carlo weight can 
become negative for three and more kinks, but in lower 
orders by f3U all diagrams are positive-definite with no 
relation to particle statistics. As a result, the system can 
be simulated at lower temperature. 

Updates. The processes chosen to update worldline 
configuration, are shown in Fig. ^. In addition, the 
time-shifting of kink should be used to speed-up cal- 
culation. For susbsystem of diagonal kinks, their cre- 
ation/annihilation and time-shifting are enough. 

Note, in case of momentum conservation p + q = r + s 
these processes can not change full momentum K of the 
system. As a result, the calculation is done in the sector 
of phase space with fixed K. This situation can be cor- 
rected by introducing fictitious kinks changing number 



ID Fermi Hubbard model. In testing purposes, we 
used the algorithm to calculate the ground state of the 
one-dimensional fermionic system with the Hamiltonian 
H^t J2 {ata^„+H.c.) + UJ2n^^n,lwitht=l,V=- 

<ij>(7 i 

1, length of the chain L — "H^ K ^ Q, number of particles 
N-^ ^ Ni = 4. The result Eq ~ -11.9523 was checked 
with exact diagonalization (£^0 = —11.952326) and real- 
space Worm algorithm. 

The average sign for both QMC algorithms as a func- 
tion of /3 is shown in Fig. |^. The graph confirms our 
assumption that sign problem is much weaker in new al- 
gorithm even for Fermi systems. With new algorithm, 
the maximal possible value of f3 in sample system is in- 
creased from 7 to 40. 

It is worthy to note that such rather good precision of 5 
digits became possible because of fixation of full momen- 
tum K = 0. The second energy level = -11.901727 
corresponding to i^T ^ is extremely close to the ground 
state, so the precision of real-space QMC is limited by 
3 digits as lowering of the temperature is prevented by 
sign problem. 

Fractional Quantum Hall Effect. In studying Wigner 
Crystallization in Fractional Quantum Hall Effect, the 
fermionic Hamiltonian of view (^) was obtained and 
then analysed using exact diagonalization method 0. 
Now the QMC algorithm described here, was applied to 
this system in order to confirm or deny the existence of 
Wigner crystal. Unfortunately, the fermionic sign prob- 
lem appeared too strong to determine the ground-state 
properties in this model. 
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FIG. 5: Average sign of MC configuration as a function of 
P = l/ksT for QMC calculations of the ID Fermi Hubbard 
model with L = 8 sites, = Ni = 4, U = -1, t = 1, 
in both real space and momentum representations. The cal- 
culations become impracticable for {Sign} < 0.01, so with 
momentum representation (and fixed full momentum K — 0) 



increases from 7 to 40. 



Bose gas with attraction in a box. The best model for 
the algorithm described in this Letter is a Bose gas in a 
box with attractive interaction. In addition to the free- 
dom of sign problem, this algorithm gets more advantage 
here, as the model can not be simulated by real-space 
QMC without errors These are caused by discretiza- 
tion of continuous real space needed to apply lattice QMC 
methods. 



II. ROTATING BOSE-EINSTEIN CONDENSATE 
WITH ATTRACTIVE INTERACTION IN ONE 
DIMENSION: HESS-FAIRBANK EFFECT AND 
MESOSCOPICS 

A. Background: macroscopic study 

Recent remarkable progress in Bose-Einstein conden- 
sation of dilute alkali gases [|| has opened up an oppor- 
tunity of studying delicate quantum phenomena in ultra- 
cold multi-atomic systems. One of intriguing set-ups is a 
system of (quasi-)one-dimensional (ID) bosons with at- 
tractive interaction - like ^Li - in a rotating ring . 
In contrast to 3D case, where the Bose-Einstcin conden- 
sate of attracting atoms becomes unstable with respect 
to a collapse above a certain threshold [which de- 
creases (vanishes) with decreasing (vanishing) the trap- 
ping potential], the ID system is unconditionally stable 
even without the trapping potential, though in a latter 
case it forms a droplet (see, e.g., Q). 



FIG. 6: Dependence of the groundstate energy E (lower 
curve) and angular momentum M (upper curve) on the rota- 
tion frequency uj. 'y — 0.1, A'^ ^ oo. 



Ueda and Leggett ^ have risen a question of whether 
the system of ID attractive bosons in a finite-size rotat- 
ing toroidal trap can remain irrotational, that is demon- 
strate the behavior typical for a superfluid - the so-called 
Hess-Fairbank (HE) effect . On the basis of their (ap- 
proximate) treatment, they arrived at the conclusion that 
the HE effect should take place in the system, provided 
the strength of interaction is below a certain threshold 
value. However, recently Berman et al. |10 have ques- 



tioned this result, arguing that the HE effect disappears 
at arbitrarily small value of attractive interaction due to 
a specific quantum instability. 

In this Section, we resolve the above-mentioned con- 
troversy by an exact numeric study of N rotating ID 
bosons with the (5-functional attractive interaction. We 
do observe the HE effect predicted by Ueda and Leggett 
(though the threshold interaction differs by a factor of 
2 from the value found in Ref. [Q). Moreover, our data 
clearly demonstrate that in the large-A^ limit the con- 
ventional methods of weakly interacting Bose gas - like 
Gross-Pitaevskii (CP) equation pl[| and Bogoliubov tech- 
nique jist - become applicable, the dimensionless small 
parameter controlling the accuracy of the approximation 
being just 1/iV. The classical-field language of CP equa- 
tion renders the issue of the presence and disappearance 
of the HE effect especially transparent. The effect per- 
sists as long as the condensate density remains uniform, 
and disappears with breaking the spatial homogeneity of 
the condensate. In a few particle system, the deviations 
from the mean- field picture are significant. In particu- 
lar, with decreasing N the HE effect gradually becomes 
indistinguishable from the generic effects of angular mo- 
mentum quantization. 

Consider N bosons of the mass m placed in the ro- 
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tating torus of radius R and cross-sectional area S = 
irr^ . At low enough temperature and with the condi- 
tion r <^ R met, the system is quasi-one-dimensional, 
and the effective ID Hamiltonian in the rotating (with 
angular frequency lo) frame reads 

^ = E (^^ - 1) ' + f E 4<^i-<,H+, ' (4) 

k k,l,q 

where the integers k,l,q stand for angular momenta, 

creates a boson with angular momentum k, n^, = cit'^k^ 
g = 2a/mRS (with a < the 3D scattering length) is the 
effective vertex of pair interaction; we use units h = 1 and 
Wc = 1, where Wc = l/2mR^ is the critical rotation fre- 
quency equal to the period of variation of the groundstate 
energy as a function of w. The groundstate of the system 
is defined by the three parameters, N, 'y = \g\{N — 1), 
and LO. 

Ueda and Leggett have analyzed the model (^) in 
the Hartree-Fock approximation in the Fock basis of an- 
gular momentum eigen states {|..., n_i, noi "-ii •••)}• They 
argued that in the limit of 7 — > 0, when no more then 
two single-particle angular momentum eigen modes sur- 
vive, their approximation is superior with respect to the 
other treatments of weakly interacting bosons. A typ- 
ical result for 7 <C 1 is shown in Fig. ^. The HF ef- 
fect - plateaus in the angular momentum curve, M(uj), - 
takes place almost at any w, except for a close vicinity of 
the critical frequency lOc- With increasing 7, the size of 
the plateaus gets smaller, and the HF effect completely 
disappears at some critical point 7 = 7c 1. At this 
point we note that 7c cannot be found accurately with 
the treatment of Ref . Q , because more than two single- 
particle angular momentum eigen modes are involved in 
the formation of the groundstate. This is immediately 
seen from the variational Hartree-Fock treatment, when 
all the particles are placed into one and the same spa- 
tially dependent single-particle state ipo (x) (non- uniform 
Bose-Einstein condensate), the wavefuction ipoix) being 
defined from the minimal energy condition, which leads 
to the GP equation 

{id/dx + cj/2f ^Po - 27r7|V'o|Vo - A^i^o - , (5) 

where fi is taken to satisfy the normalization condition 
/ |'0o(a;)P dx — 1. Solving Eq. (||) reveals the nature of 
the HF effect, as well as the mechanism of its disappear- 
ance. At w < w*(7), with w*(7) satisfying 

(w, - 2kf = 1 - 27 , (6) 

the solution ipoi^) is uniform and thus irrotational. At 
uj > LO^('y), the rotational symmetry of the problem 
breaks down: The density becomes non-uniform, and the 
rotation of the density profile gives rise to the increase of 
the angular momentum with growings. At 7 > 7c = 1/2 
the minimal-energy solution is non-uniform even without 
the rotation and the HF effect totally disappears. 
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FIG. 7: Function M{lo) at 7 = 0.45, as predicted by GP 
equation (upper solid curve). Lower solid curve represents 
LJ, (7), according to Eq. (^. Dotted curve is the prediction of 
Ref. |] for Lo*{j). 

A simple and instructive way of arriving at the relation 
(^) is the Bogoliubov treatment jl^ for the elementary 
excitation spectrum in assumption of the uniform con- 
densate. Considering the energy, ei, of the first excited 
state, one observes that at w > w*, the uniform conden- 
sate is thermodynamically unstable: ei < 0. At 7 > 1/2 
the energy ei becomes imaginary, indicating dynamical 
instability of the homogeneous solution at any lu. 

At 7 ^ 1 the solution of Eq. (H) is strongly non- 
uniform and corresponds to a rotating condensate droplet 
- bright soliton 0, 

To obtain the GP prediction for the angular momen- 
tum as a function of we solved Eq. (|) by numeri- 
cally minimizing GP energy functional, see Figs. 0, ||. 
At 7 ^ 1, the GP results coincide, up to higher-order 
corrections, with those of Ref. At 7 ~ 1, however, 
deviations become significant. 



B. Simulation for finite A'^ 

As we have already mentioned, the very existence of 
the HF effect (and thus the applicability of the weakly 
interacting gas treatments) has been questioned recently 
by Herman et al. [0. This conclusion of Ref . ^ seems 
rather strange and counter-intuitive in view of the known 
exact results for the ID attractive bosons in the non- 
restricted geometry. The small parameter that guaran- 
tees applicability of the mean-field approach is just the 
inverse number of particles [ p^ . A priori we do not see 
how the finite system size can qualitatively change the 
situation. 

In view of this controversy, as well as keeping in mind 
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FIG. 8: Gross-Pitaevskii equation results for the function 
M{ijj) at different 7's. 



the fact that mesoscopic results are interesting in their 
own value, we performed an exact numeric study of the 
model (|4|) for different numbers of bosons. 

Analytic summation of diagonal kinks contribution due 
to point interaction U^qrs = C^^p+g.r+s increases full 
energy by constant Ai? = Uo{2N'^ — N), and parti- 
cle energy obtains term depending on occupation ik = 
{k-Lu/2) -UoUk- 

After typical update, the weight of MC configuration 
acquires the factor ~ g-5S(r-ro) ]y[ain difference from 
lattice models is that SE can be unlimitedly large, so 
with large SE, the weight of new configuration would be 
too small to appear in simulation and the efficiency of 
the update approaches zero. Therefore we must choose 
for update different places with different probability to 
make simulation effective enough. 

Note, for old — {k~uj/2)^, there exist simple re- 
lation 5E = 2q{q + A) with A characterizing the place 
of update. Therefore we can neglect second term —UoUk 
in new particle energy e, and choose q using simple an- 
alytics, given in Appendix ^ for all types of updates. 
Though for small fc, both terms in new particle energy 
ik are comparable (for in this model Uq = 7/2 (iV — 1), 
7 < 1), this trick is aimed mostly for resonable choos- 
ing of large momentum k where (nj.) ~ making second 
term rather virtual. 



C. Results 

We traced the evolution of the system properties with 
N variying from 2 to 100, at different w's and 7's. 
The simplest characteristics that we studied was the 
groundstate energy, which we calculated by simulating 
the groundstates in different angular momentum sectors. 
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FIG. 9: Groundstate energies in different angular momentum 
sectors for Ai' = 10 at 7 = 1/2. The bars are the Monte 
Carlo data with errors. The solid curves are the parabolic 
fits. Integers stand for the values of the angular momenta. 



with subsequent selecting the global minimum, see Fig. |^. 
This procedure yields also the curve M{uj). For a finite- 
system this curve is essentially stepwise due to the quan- 
tization of angular momentum. Comparing this curve at 
large enough N to the Gross-Pitaevskii solution, we find 
an excellent agreement, see Fig. 

To quantitatively trace the difference between GP and 
Monte Carlo results, we compared corresponding answers 
for the groundstate energies at different A^'s. We found 
that starting from N Ri 10 the deviation between the two 
results scales as 1/A^, which confirms that in the A^ 
00 limit the Gross-Pitaevskii equation yields a perfect 
description of the groundstate properties. 

In Fig. ^ we present the M{uj) curves for 7 = 0.2 < 
7c — 1/2, at different particle numbers. Once again we 
see an excellent agreement with the GP equation at large 
A^. The HF plateau is unambiguously revealed. Note, 
that at TV < 5 the HF plateau at zero momentum is indis- 
tinguishable (by its size) from the rest of the quantized- 
momentum plateaus. 

To get an insight into the inner structure of the ground 
state, we calculate the two-particle densityn correlator 
K{x) = ^+{x'^+{x' + x)-^{x' + x)^{x'))^, /N{N-1). 
In Figs, naand Owe present K{x) {uj = 0) for 71 — 0.25 
and 0.75, that is for uniform and non- uniform (in the 
macroscopic limits) cases, respectively. At small enough 
A^ there is no qualitative difference between the two cases. 
At large A^ the difference is clearly seen. Once again note 
an excellent agreement with the GP equation. 

Summarizing, we developed a novel Quantum Monte 
Carlo algorithm based on the generic principles of the 
Diagrammatic Monte Carlo approach Q . The algorithm 
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FIG. 10: Groundstate angular momentum M as a function of 
the rotation frequency w for AT = 10 at 7 = 1/2 (solid line). 

Absolute error is on the order of 10~^. Dotted curve is the 
N ^ 00 limit obtained by solving Gross-Pitaevskii equation. 



FIG. 12: Density-density correlator K{x) = 0) at 7 = 
0.25 < 7c = 1/2 for iV = 80, 10,4,2. The correlations grow 
up with decreasing N. 
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FIG. 11: Groundstate angular momentum M as a function of 
rotation frequency cu for N = 5, 10, 20 at 7 = 0.2 < 7c = 1/2. 
Absolute error is on the order of 10~^. The dashed line is the 
result of the Gross-Pitaevskii equation. 



samples exact diagrammatic expansion (in terms of the 
pair interaction) of the imaginary-time evolution opera- 
tor in the m,omentum representation. In comparison with 
usual real-space methods, this algorithm is more efficient 
in simulation of weakly interacting systems. Moreover, 
it can be apphcd to the models which can not be stud- 
ied by real-space methods. The sign problem is weaker 
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FIG. 13: Density-density correlator K{x) (w = 0) at 7 = 
0.75 > 7c = 1/2 for N = 80, 10, 4, 2. Sohd curve corresponds 
to the Gross-Pitaevskii equation. The particle positions are 
correlated even in the macroscopic limit indicating the non- 
uniformity of the groundstate. 



here making possible to reach temperatures low enough 
to study ground state. In the case of attractive pair po- 
tential, all diagrams are positive- definite and the method 
is very efficient, allowing to simulate up to 100 and more 
particles. 

Using exact Quantum Monte Carlo algorithm, the 
groundstate properties of one-dimensional bosons with 

attractive (5-functional interaction in a rotating toroidal 
trap are studied. The Hess-Fairbank effect - absence of 
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a response to the trap rotation - is observed in a cer- 
tain area of the parameter space. The fact that in the 
N ^ oo Umit the Gross-Pitaevskii equation yields a per- 
fect description of the groundstate properties is proved. 

The author is grateful to Profs. V.A. Kashurnikov, 
B.V. Svistunov, and N.V. Prokof'cv for drawing my at- 
tention to this problem and numerous helpful discussions. 
This work was supported by the Russian Foundation for 
Basic Research. 
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APPENDIX A: WEIGHTING THE TIME OF 
NEW KINK 

With adding new kink in time r in range 
(Tmin ■ ■ ■ Tmax) , the Statistical weight of MC configura- 
tion receives a factor ^ g-Q(T-r„,i„)^ where Q denotes 
change of worldline energy in the range (r,„i„ . . . t) after 
update. 

Time r is determined with probability 



using relation 



P(Q,r, 



max ' mm ■) -^^/i 



R), 



where R is the random number uniformly distributed be- 
tween and 1. 

For convenience, we define the following functions: 



ZiiQ,T) 



usually, 
if IQI ^0, 
if Q -l-oo, 
if Q ^ —CO, 



(Al) 



p(g,T,i?) 



(-Ml±M^^^ usually, 
RT, if \Q\ -> 0, 

if Q — > +00, 
if Q — *■ — oo. 



ilni?, 
T-^lni?, 



(A2) 

Limiting cases should be realized separately to avoid al- 
gorithm inefficiency and precision loss. 



APPENDIX B: UPDATES IN DETAILS 

a. Creation/ Annihilation of Two Kinks 

This process is shown in Fig. 
Two kinks are created as follows: 

1 . Time of first kink tq from to /3 is chosen in random 
way with probability 

2. Two not-empty world lines are found going through 
To with momenta ki and ^2, including possible case 
ki = /c2. Let us denote the probability to choose 
these world lines W{ki, k2). 



FIG. 14: Creation/ Annihilaton of Two Kinks. 



3. The second kink will be created at the time ti de- 
termined using Appendix 

Tl = To + P(2gf(gf + A), Tmaa; - Tq). 

With increasing fca and k^^ their average occupa- 
tion becomes exponentially small, so Tmax remains 
intact. This fact allows us to take into account all 
possible kz,ki in single process. The analytics for 
choosing q is given in Appendix 

4. The probability to accept this process is found from 
the balance equation (Bl). 



The scheme of annihilation: 

1. A random kink having parameters tq, /ci, fc2, 
5, Ueff, neff is taken with probability -^^ — = 

1 

Qald+2- 

2. The kink chosen and its right neighbour (the direc- 
tion is fixed to remove polyvalence) must form the 
"kink-antikink" pair. 

3. The probability to choose q in the direct process is 

W{q) 

4. The annihilation probability pj^'^'^^ of this pair is 
determined from the balance equation (Bl). 



WoldW{ki,k2) 



W{q) At Are^- 

Z3 ~p 



-P 



(acc) 



(Bl) 



Woidi-ATUeSfn.ffYe 



Q + 1 



(acc) 



(Pay attention to the possibility of pair annihilation in 
two different ways when all four momenta ki, ...,ki are 
untouched by other kinks, see Fig. |l5|. To remove this 
duality we disable the annihilaton when ti < tq). 



b. Kink Expanding/Contracting - particle version 
The direct process (see Fig. ^6|) is done as follows: 
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(b) 



(c) 
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FIG. 15: About possibility to annihilate a pair of kinks in two 
ways resulting in different configurations. 



ks 
•ks 



() A 



^0 ^1 



FIG. 16: Process of Kink Expanding/Contracting - particle 
version. 



1. A random kink having parameters Tq, fci, fc4, Ui, 
ni is chosen with probabiHty ^ . 

2. Choice of q and ti is done similarly to the cre- 
ation of two kinks with probabilities W{q)/Z3 and 
exp(— Ai?(Ti — tq))/Zi, respectively. 

3. The probability to accept this update is determined 
using balance equation (B2). 



Inverse transformation: 



1. A random kink having parameters tq, fci,...,fc4, 
U2, Us, n2, ns is chosen with probability — = 



2. The kink chosen must, with its right neighbour, 
form the pair which can be contracted into single 
kink. 

3. The probability to choose this q in direct process 
equals W{q)/Z3. 

4. The probability to to accept this update is deter- 
mined from the balance equation (B2). 



FIG. 17: Kink Expanding/Contracting - "hole" version. 

c. Kink Expanding/Contracting - "hole" version 

In contrast to previous update, this process needs mo- 
menta /cs, ke to be occupied. However, in each moment 
To only limited set of occupied pairs (k^jke) exists with 
^5 + ^6 = ^1 + ^2 (see Fig. |l^), i.e. fitting this update. 
Therefore the momenta should not be weighted using for- 
mulas of Appendix Otherwise, we choose a random 
pair from rather small set. 

The direct update is made as follows: 

1. A random kink having parameters tq, fci, fc4, Ui, 
ni is chosen with probability ^. 

2. All pairs {k^, ke) fitting this process are determined; 
for each pair, the parameters T^ax, ?^2, "-s, U2, U3, 
AE are found. 

3. Each pair (k^^ke) is chosen with probability 
W{k5,ke) ~ U2U3n2n3Zi{AE 



4. The balance equation ( JBSj ) is used to determine the 
probability to accept this update. 

The inverse process is made as follows: 

1. A random kink having parameters tq, fci, k2, k^, fcgj 
U2, U3, n2, is chosen with probability — = 

1 

2. The kink chosen and its right neighbour must form 
the pair which can be contracted into single kink 
fitting the direct "hole" update. 

3. The probability 1/^(^5, fcg) to choose this pair 
(fcs, fcg) in direct process, is determined. 

4. The probability to accept this update is determined 
from the balance equation (B3). 



1 W(a) Are^---) / ^ 
W^~,„4-ATf/ini)-^— = (B2) 



common 

Wcommon{-ATyU2U3n2n3e 



(■■■). 



1 



(0 + 1) 



[ace] 



w, 



common 

w, 



{^ATUin,)^W{k5,ke)^pt'^^ = (B3) 



^ common 



{~ATfU2U3n2n3e 



(■■■). 



1 



(0 + 1) 



(acc) 
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FIG. 18: Entangling of a Kink with a Worldline. In the 
case shown here, the upper parts of left and right kinks are 
touched. 
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FIG. 20: To the process of Shifting the Kink Through An- 
other. In simulation of periodical chain, two kinks can occupy 
all the same momenta though not forming the "kink-antikink" 
pair. L denotes the length of a chain. 



tical. The scheme is as follows: 

1. A random kink having parameters tq, fci,...,A:4 is 
chosen with probabiUty ^; 



FIG. 19: Shifting the Kink Through Another. Kinks exchange 
their times; the occupation numbers are changed for world- 
lines with momenta common for both kinks. 



2. The nearest kink on the right with at least one the 
same momentum, is found; 

3. The parameters AEbefore, ^Eafter, ribefore, riafter 

are determined; 



d. Entangling of a Kink with a Worldline 

This process can not be made of the processes de- 
scribed above. More visually, this is the only process 
able to change the winding numbers of worldlines. There 
exist many versions of this update. The case touching 
upper parts of left and right kinks is shown in Fig. 

The momentum fcg should be occupied at tq, thus there 
is only limited set of pairs {kc,^kQ) fitting the process of 
entangling. Therefore, the schemes of direct and inverse 
updates are similar to that for the previous transforma- 
tion. The only correction should be made is to choose 
upper or lower parts of kinks, if present. 



e. Shifting the Kink Through Another 

Carrying out shift in time, one must take into account 
possible collisions with neighbouring kinks. While simple 
shift is enough for simulating Fermi systems, Bose case 
should incorporate shift "through" neughbouring kink 
(Fig. |l^). In this process, the kinks are exchanged by 
time: t^^'"'^ = rf'^, t^^'" = t^'''. 

The number K of common momenta can be 1 to 4. The 
case = 4 is usually associated with the "kink-antikink" 
pair, so the shifting through the kink can be replaced by 
the combination of Annihilation and Creation. However, 
another case shown in Fig. ( pO[ ) can appear in simulation 
of Bose chain. 

Due to symmetry, direct and inverse processes are iden- 



4. The probability to accept t he u pdate is determined 
from the balance equation (B4). 



or, finally. 



1 



A-Ebe/ore (t2 -Ti ) _|_ (acc) 



(B4) 



common 



a 



{ace) 



(acc) 



ngfter ^-(A£„ffer-A£b,f„,e)(r2-Ti) 
'^before 



APPENDIX C: CHOOSING MOMENTUM 
ANALYTICALLY 

We use the fact 5E — 2q{q + A), where A is fixed in 
choosing the place of action, and q > —A/2. The value 
of Zi{SE{q),T) as a function of q, is shown in Fig. |2|. 
It can be approximated by following piecewise function 



(ZiiEiq),T), ifq<0 
Zi^ W{q) = < const ^T, if < q < q* , 



(CI) 



making possible to determine q analytically. 
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W(q) 



correspond to three pieces of approximation (CI). Values 



Z^-\ define boundary values = 4r ^^d 

R(+) = 1 - 




To choose the random number R uniformly dis- 
tributed between and 1. is taken. Initially, we de- 
termine the range {q < Q < q < q* , q > q*), then q 
according to the range found. When R < R^-\ we find 

q'<q 

the value q < 0, making J2 W{q') greater than RZ^. 



Other cases allow analytical relation: q = 1+q* 
if < i? < and 



-A/2 



10 15 20 



FIG. 21: Zi{T,2q{q + A)) as a function of q (solid curve) j ^, A = 0, 

and its piecewise approximation W(q) (dotted curve) used to '^^1 A>0 

choose q analytically. I -1 ' ' 



The choice of q* is based on the relation 1/E{q*) sa T: 



if A = 0, 



where r = -ff|n^, if i? > R'^' 



(C2) 



The possibility to choose given q equals where 
Z3 = Z(-) + Z(0) and 



Z(-)= ^ W(g), 
9* 



g>0 



C30 ^ 

~ V - 

2g(g + A) 

q>q 



if A = 0, 



^Infl + A] ifA>0 



(C3) 
(C4) 



(C5) 



The scheme modifications for cases A > — 1 (without 
range q < 0) and q* — (without range < q < q*) are 
straightforward. 

Note. In summation Z'^', the care should be taken 
in managing overflow. The exponential index should not 
be limited by any value. Otherwise the relations between 
values of Zi would be broken, making most probable mo- 
menta q < with nonequal energies {q — uj/2(jjc)'^ to have 
similar weights. 



Some way to avoid this, is to decrease the energy dif- 
ference SE = - E"'-'^ = 2q{q + A) by some value 
into safe region. 
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